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Abstract: Hidden Markov Models (HMMs) are widely used probabilistic models, 
particularly for annotating sequential data with an underlying hidden structure. Patterns 
in the annotation are often more relevant to study than the hidden structure itself. A typical 
HMM analysis consists of annotating the observed data using a decoding algorithm and 
analyzing the annotation to study patterns of interest. For example, given an HMM modeling 
genes in DNA sequences, the focus is on occurrences of genes in the annotation. In this 
paper, we define a pattern through a regular expression and present a restriction of three 
classical algorithms to take the number of occurrences of the pattern in the hidden sequence 
into account. We present a new algorithm to compute the distribution of the number of 
pattern occurrences, and we extend the two most widely used existing decoding algorithms 
to employ information from this distribution. We show experimentally that the expectation 
of the distribution of the number of pattern occurrences gives a highly accurate estimate, 
while the typical procedure can be biased in the sense that the identified number of pattern 
occurrences does not correspond to the true number. We furthermore show that using this 
distribution in the decoding algorithms improves the predictive power of the model. 
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1. Introduction 

A Hidden Markov Model (HMM) is a probabilistic model for sequential data with an underlying 
hidden structure. Because of their computational and analytical tractability, they are widely 
used especially in speech recognition [1-3], image processing [4] and in several applications 
in bioinformatics; e.g., modeling of proteins [5-7], sequence alignment [8-10], phylogenetic 
analysis [11,12] and identification of coding regions in genomes [13,14]. 

Patterns in the hidden structure are, however, often more relevant to study than the full hidden 
structure itself. When modeling proteins, one might be interested in neighboring secondary structures 
that differ, while for sequence alignments, the pattern could capture specific characteristics, such as long 
indels. In phylogenetic analysis, changes in the tree along the sequence are most relevant, while when 
investigating coding regions of DNA data, patterns corresponding to genes are the main focus. 

Counting the number of occurrences of such patterns can be approached (as in the methods based 
on [15]) by making inferences from the prediction of a decoding algorithm; e.g., the Viterbi algorithm 
or the posterior- Viterbi algorithm. As we show in our experiments, this can give consistently biased 
estimates for the number of pattern occurrences. A more realistic method is presented in [16], where the 
distribution of the number of pattern occurrences is computed by means of Markov chain embedding. 
To our knowledge, this is the only study of patterns in the hidden sequence of HMMs. The problem 
of pattern finding in random sequences generated by simple models, such as Markov chains, has been 
intensively studied using the embedding technique [17-19]. 

We present a fundamentally different approach to compute the distribution of the number of pattern 
occurrences and show how it can be used to improve the prediction of the hidden structure. We 
use regular expressions as patterns and employ their deterministic finite automata to keep track of 
occurrences. The use of automata to describe occurrences of patterns in Markov sequences has been 
described previously in [18,20,21]. However, analyzing pattern occurrences in the hidden structure of 
HMMs by means of automata has not been done before. We introduce a new version of the forward 
algorithm, the restricted forward algorithm, which computes the likelihood of the data under the hidden 
sequence containing a specific number of pattern occurrences. This algorithm can be used to compute the 
occurrence number distribution. We furthermore introduce new versions of the two most widely used 
decoding algorithms, the Viterbi algorithm and the posterior- Viterbi algorithm, where the prediction 
is restricted to containing a certain number of occurrences, e.g., the expectation obtained from the 
distribution. We have implemented and tested the new algorithms and performed experiments that show 
that this approach improves both the estimate of the number of pattern occurrences and the prediction of 
the hidden structure. 

The remainder of this paper is organized as follows: we start by introducing Hidden 
Markov Models and automata; we continue by presenting our restricted algorithms, which we then 
validate experimentally. 
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2. Methods 



2.1. Hidden Markov Models 



A Hidden Markov Model (HMM) [3] is a joint probability distribution over an observed sequence 
Vi-.T = yiy2---yT ^ O* and a hidden sequence Xi.t = XiX2--.xt G V.*, where O and V. are 
finite alphabets of observables and hidden states, respectively. The hidden sequence is a realization 
of a Markov process that explains the hidden properties of the observed data. We can formally 
define an HMM as consisting of a finite alphabet of hidden states T-L = {hi, h2, . . . , h]\f}, a finite 
alphabet of observables O = {01,02, . . . ,om}, a vector 11 = {nhJi<i<N, where tt/j. = P(Xi = hi) 
is the probability of the hidden sequence starting in state hi, a matrix A = {ahi^hj}i<i.j<N, where 
0'hi,hj = P(-^t = hj I Xt-i = hi) is the probability of a transition from state hi to state hj, and a matrix 
B = {h„o,}i<i<N , where bh^^o^ = '^iYt = Oj \ Xt = hi) is the probability of state hi emitting oj. 

Figure 1 shows an HMM designed for gene finding. The observed sequence is a DNA sequence over 
the alphabet O = {A, C, G, T}, while the hidden states encode if a position is in a non-coding area (A^) 
or is part of a gene on the positive strand (C) or on the negative strand (R). The model incorporates the 
fact that nucleotides come in multiples of three within genes, where each nucleotide triplet codes for an 
amino acid. The set of hidden states is "H = {A^} U {Ci, Ri\l < i < 3}. In practice, models used for 
gene finding are much more complex, but this model captures the essential aspects of a gene finder. 

Figure 1. A Hidden Markov Model (HMM) for gene prediction. Each box represents 
a hidden state, and the numbers inside are the emission probabilities of each nucleotide. 
Numbers on arcs are transition probabilities between hidden states. 



C3 




C2 




Ci 




N 




R3 




R2 




Ri 


A: 0.30 


1 


A: 0.20 


1 


A: 0.40 


0.05 


A: 0.25 


0.05 


A: 0.20 


1 


A: 0.30 


1 


A: 0.15 


C: 0.25 




C: 0.35 




C: 0.15 




C: 0.25 




C: 0.40 




C: 0.20 




C: 0.30 


< 


< 


< 


> 


> 


> 


G: 0.25 




G: 0.15 




G: 0.20 




G: 0.25 




G: 0.30 




G: 0.30 




G: 0.20 


T: 0.20 




T: 0.30 




T: 0.25 




T: 0.25 




T: 0.10 




T: 0.20 




T: 0.35 



0.9 



0.9 



HMMs can be used to generate sequences of observables, but their main application is for analyzing 
an observed sequence, yi-^. The likelihood of a given observed sequence can be computed using the 
forward algorithm [3], while the Viterbi algorithm [3] and the posterior- Viterbi algorithm [22] are used 
for predicting a corresponding hidden sequence. All these algorithms run in O (TN'^) time, using 
O (TN) space. 

2.1.1. The Forward Algorithm 

The forward algorithm [3] finds the probability of observing yi,T by summing the joint probability of 
the observed and hidden sequences for all possible sequences, xi:t- This is given by: 

T 

P(|/1:T, Xi-.t) = T^xJ)x^,yr JJ a^.t^^,xtbxt,yt (1) 

t=2 
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P(2/1:t) = 5ZP(|/1:T,Xi;t) 



(2) 



where Equation (1) is the mukiplication of the probabilities of transitions and emissions, which explain 
observing yi,T with xi,t as the hidden sequence: the HMM starts in state xi and emits yi from xi, and 
for alH = 2, . . . , T, it makes a transition from state xt^i to Xt and emits yt from xt- 

The forward algorithm finds Equation (2) by recursively filling up a table, a, with values 
at{xt) = F{yi:t,xt) = Y.x,..t-i P(2/i:t' ^i:t) being the probability of observing yi-t and being in state 
Xt at time t. The recursion is: 



and, finally, P(i/i:r) = Y.i c^rihi). 

2.1.2. The Viterbi Algorithm 

The Viterbi algorithm [3] finds the sequence of hidden states, xi:t, that maximizes the joint probability 
of the observed and hidden sequences (1). It uses the same type of approach as the forward algorithm: a 
new table, uj, is defined by u)t{xt) = max^.^ {T{yi,t, Xi;^)}, the probability of a most likely decoding 
ending in xt at time t, having observed yi.t- This can be obtained as follows: 



After computing w, a most likely sequence of hidden states is retrieved by backtracking through the 
table, starting in entry argmax^. {a;r(/ij)}. 

2.1.3. The Posterior Decoding Algorithm 

The posterior decoding [3] of an observed sequence is an alternative to the prediction given by the 
Viterbi algorithm. While the Viterbi algorithm computes the decoding with the highest joint probability, 
the posterior decoding computes a sequence of hidden states, xi-t, such that Xt = argmax^. {7* (/ij)} has 
the highest posterior probability 7t(/ii) = P(/?.i | yi-.r)- If we let l3t{h,i) = F{yt+i;T \ hi), we have: 



ai{hi) 




(3) 



uJt{hi) = bh„yt ■ max {ut-i{hj) at^^h,} 



(4) 



F{h„y^.,t)F{yt+,.,T \ K) 
P(yi:T) 



(5) 



The backward algorithm [3] gives (3t{hi) by using a recursion similar to Equation (3). Thus, to 
compute the posterior decoding, we first fill out at{hi) and (3t{hi) for all t and i and then compute the 
decoding by Xt = argmax^^ {7t(^i)}- 
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2.1.4. The Posterior- Viterbi Algorithm 

The posterior decoding algorithm often computes decodings that are very accurate locally, but it may 
return syntactically incorrect decodings; i.e., decodings with transitions that have a probability of zero. 
The posterior- Viterbi algorithm [22] corrects for this by computing a syntactically correct decoding 
x\.rp = argmax^^^g^^ |nLi7i(^*)} "^^^^ '^he highest posterior probability, where Ap is 
the set of syntactically correct decodings. To compute this, a new table, 7, is defined by 
Itixt) = max^j {nLi lt{xt)}, the maximum posterior probability of a decoding from Ap ending 

in Xt at time t. The table is filled using the recursion: 

lt{hi) = -ftihi) ■ max {lt-i{hj)} 

After computing 7, a decoding from Ap with the highest posterior probability is retrieved by 
backtracking through 7 from entry argmax^^ {irihi)}. We note that, provided that the posterior 
decoding algorithm returns a decoding from Ap, the posterior- Viterbi algorithm will return the same 
decoding. 

2.2. Automata 

In this paper, we are interested in patterns over the hidden alphabet 7i = {hi, h2, . . . , h^} of an 
HMM. Let r be a regular expression over "H, and let FA-}{(r) = (Q, "H, go, ^, [23] be the deterministic 
finite automaton (DFA) that recognizes the language described by (/ii | /i2 | . . . | /iAr)*(r), where Q is 
the finite set of states, go ^ Q is the initial state, A C Q is the set of accepting states and 5 : Q xT-L ^ Q 
is the transition function. Fy4^(r) accepts any string that has r as a suffix. We construct the DFA 
EA'^{r) = {Q, v., go, A, 6e) as an extension of Fyl^(r), where 6e is defined by: 

WqeQ\A, When, 6E{q,h) = 6{q,h) 
WqeA, When, 6E{q,h) = 6{qo,h) 

Essentially, EAy^{r) restarts every time it reaches an accepting state. Figure 2 shows F^-^(r) and 
EAn{r) for r = (NCi) \ (RiN) with the hidden alphabet n = {N} U {Ci,Ri | 1 < i < 3} of the 
HMM from Figure 2. Both automata have Q = {0, 1, 2, 3, 4}, go = 0, and A = {2, 4}. State 1 marks the 
beginning of NCi, while state 3 corresponds to the beginning of RiN. State 2 accepts NCi, and state 4 
accepts -RiA^. As C2, C3, R2 and R3 are not part of r, using them, the automaton restarts by transitioning 
to state 0 from all states. We left these transitions out of the figure for clarity. The main difference 
between Fy4-^(r) and EAy^(r) is that they correspond to overlapping and non-overlapping occurrences, 
respectively. For example, for the input string, RiNCi, FAy^ (r) first finds RiN using state 4, from which 
it transitions to state 2 and matches NCi. However, after EAy^(r) recognizes -RiA^, it transitions back to 
state 0, not matching NCi. The algorithms we provide are independent of which of the two automata is 
used, and therefore, all that remains is to switch between them when needed. In our implementation, we 
used an automata library for Java [24] to obtain FA'y{r), which we then converted to EA'y{r). 
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Figure 2. Two automata, F^^(r) and EA^ir), for the pattern r = (NCi) \ {RiN), 
n = {N} U{Ci,Ri \ 1 < i < 3}, Q = {0,1,2,3,4}, go = 0, and A = {2,4}. States 
1, 2 and 3, 4 are used for matching sequences ending with NCi and RiN, respectively, as 
marked with gray boxes. The two automata differ only with respect to transitions from 
accepting states: the dotted transition belongs to Fyl^(r) and the dashed one to EA-^^r). For 
clarity, the figure lacks transitions going from all states to state 0 using C2, C3, R2 and R3. 




3. Results and Discussion 

Consider an HMM as defined previously, and let r be a regular expression over the alphabet of 
hidden states, "H. We present a modified version of the forward algorithm to compute the distribution 
of the number of occurrences of r, which we then use in an adaptation of the Viterbi algorithm and the 
posterior- Viterbi algorithm to obtain an improved prediction. 



3. 1. The Restricted Forward Algorithm 

Let Or{xi:T) be the number of matches of r in xi-t- We wish to estimate Or by using its probability 
distribution. We do this by running the HMM and Fyl^(r) in parallel. Let FAfi{r)t be the state in 
which the automaton is after t transitions, and define at(xt,k,q) = '^[xi.t-i-Or{xvt)=k}^iy'i^-t^ ^^-t: 
FA-^{r)t — q) to be the entries of a new table, a, where k = 0,...,m and m < T is the 
maximum number of pattern occurrences in a hidden sequence of length T. The table entries are the 
probabilities of having observed yi,t, being in hidden state Xt and automaton state q at time t and having 
seen k occurrences of the pattern, corresponding to having visited accepting states k times. Letting 
5~^{q,hi) = {q' \ 6{q',hi) = q} be the automaton states from which a transition to q exists, using 
hidden state hi and 1 being the indicator function, mapping a Boolean expression to one if it is satisfied 
and to zero otherwise, we have that: 

at{xt, k,q)^ P(yi:t, Xi.,t, FAn{r)t = q) (8) 

Xl:t-1 ■■ 

Or{xv.t-i)=k-l{qeA) 

P(yi:t, xi..t, FAH{r)t = q)= J] P(|/i:t, xi..t, FAnir)t-i = q') 

= hxt,yt(^xt-i,xt Y xi.,t-i, FA'n{r)t-i = q') 

g'e5-'^{q,xt) 
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Using Equations (8) and (9), the recursion for a becomes: 

{t{q^A) ifA; = 0 

l{qeA) iik = l 

0 otherwise (10) 

at{hi,k,q) = bh^^y,^ah^^h, ^ at-i{hj,k - t{q e A),q') 

These probabilities now allow for the evaluation of the distribution of the number of occurrences, 
conditioned on the observed data: 

^ , , P(fc occurrences of r, wi-t) 
rik occurrences of r ui-t) = ;^;7 r ' — 

from which the expectation can be computed. The likelihood of the data can be obtained from either the 
forward or restricted forward algorithm, F{yi,T) = Xlj ct(^i) = ^ikq '^rihi, k, q). 

The a table has T ■ N ■ m ■ \Q\ entries, and the computation of each requires O {N\Q\), leading to 
a O {TN'^m\Q\'^) running time and a space consumption of O {TNm\Q\). In practice, both time and 
space consumption can be reduced. The restricted forward algorithm can be run for values of k that are 
gradually increasing up to kmax for which P(at most kmax occurrences of r | yi-.r) is greater than, e.g., 
99.99%. This k^ax is generally significantly less than m, while the expectation of the number of matches 
of r can be reliably calculated from this truncated distribution. The space consumption can be reduced 
io O {N\Q\), because the calculation at time t for a specific value, k, depends only on the results at time 
t — 1 for A; and k — 1. 



3.2. Restricted Decoding Algorithms 

The aim of the restricted decoding algorithms is to obtain a sequence of hidden states, xi.t, for which 
Or{xi,T) G [/, u], where / and u are set to, for example, the expected number of occurrences, which can 
be calculated from the distribution. The restricted decoding algorithms are built in the same way as the 
restricted forward was obtained: a new table is defined, which is filled using a simple recursion. The 
evaluation of the table is followed by backtracking to obtain a sequence of hidden states, which contains 
between / and u occurrences of the pattern. The two restricted decoding algorithms use O {TNu\Q\) 
space and O {TN'^u\Q\^) time. 

The entries in the table for the restricted Viterbi algorithm contain the probability of a most likely 
decoding containing k pattern occurrences, ending in state Xt and automaton state q at time t and having 
observed yi.t, uJt{xt, k, q) = max{^.j^^_^ ;o,(xi.t)=fc} {^{Vi-.t, Xi.t, FA-H{r)t = q)} with k = 0,...,u. The 
corresponding recursion is: 

{l{q^A) ifk = 0 
l{q e A) ifA; = l 
0 otherwise (12) 

ujt{hi, k, q) = bh„yt ■ max <^ ah,,h, ■ max {ut-i{hj, k-l{q e A), q')} \ 

j L g'G(5-i((3,/ii) I 
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and the backtracking starts in entry argmax-^gj;^] ^ {c2;r(/ij, k, q)}. 

For the restricted posterior- Viterbi algorithm, we compute the highest posterior probability of a 
decoding from Ap containing k pattern occurrences, ending in state Xt and automaton state q at time 
t and having observed yi.t, ^t{xt, k, q) = max{^^.^_^,o,{xi:t)=k} {'^ixi-.t, FAy,{r)t = q \ yi-.r)}- We have: 



t{q^A) ifk = 0 

^i{hi,k,q)=-fi{h,)l{qoe6-\q,h,))- { l{q E A) if /c = 1 

0 otherwise 



(13) 



{q'e5-^{q,hi) ) 

The backtracking starts in entry argmax-^gj^^j^ {irihi, k, q)}. 



^t{hi, k,q) = ^t{hi) ■ max 

{hj:ah.h.>0} 



3.3. Experimental Results on Simulated Data 

We implemented the algorithms in Java, validated and evaluated their performance experimentally as 
follows: We first generated a test set consisting of 500 pairs of observed and hidden sequences for each 
length L = 500, 525, . . . , 1,500 from the gene finder in Figure 2. As the HMM is used for finding genes, 
an obvious choice of pattern is r = (A^Ci) | (RiN), corresponding to the start of a gene. For each of 
the sequences, we estimated the number of overlapping pattern occurrences with the expected number 
of pattern occurrences, computed using the restricted forward algorithm, which we then used to run the 
restricted decoding algorithms. We also computed the prediction given by the two unrestricted decoding 
algorithms for comparison. 



Figure 3. Normalized difference. 



1, between the true number of pattern 



true value 

occurrences, the number given by the two unrestricted decoding algorithms and the expected 
number of pattern occurrences computed using the restricted forward algorithm. For each 
sequence length, we show the median of the normalized differences in the 500 experiments, 
together with the 0.025 and 0.975 quantiles, given as error bars. 
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3.3.1. Counting Pattern Occurrences 

Figure 3 shows the power of the restricted forward algorithm and the two unrestricted decoding 
algorithms to recover the true number of pattern occurrences. For each given length, we computed 
the normalized difference, -^^iimste — ^j^^j plotted the median over the 500 sequences together with 
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the 0.025 and 0.975 quantiles, given as error bars. As Figure 3 shows, the expectation gives a very 
good estimate. The decoding algorithms' performances are significantly lower, and they always give 
underestimates. This may be partly due to the model structure, where transitions to and from coding 
regions have low probability, leading to few, but long, genes. 

3.3.2. Quality of Predictions 

We compared the predictive power of the two decoding algorithms in the original and restricted 
versions, using the expectation for the number of pattern occurrences. For each length in the test set, 
we measured the quality of each method, both at the nucleotide and gene level, following the analysis 
in [25]. Because the measures we use require binary data, we first converted the true hidden sequence 
and decoding to binary data containing non-coding areas and genes, by first considering the genes on the 
reverse strand as non-coding and calculating the measures for the genes on the direct strand, and vice 
versa. The final plotted measures are the averages obtained from the two conversions. 

Nucleotide level. To investigate the quality at the nucleotide level, we compared the decoding and the 
true hidden state position by position. Each position can be classified as a true positive (predicted as part 
of a gene when it was part of a gene), true negative (predicted as non-coding when it was non-coding), 
false positive (predicted as part of a gene when it was non-coding) and false negative (predicted as 
non-coding when it was part of a gene). Using the total number of true positives {tp), true negatives 
{tn), false positives (Jp) and false negatives (/n), we calculated the sensitivity, specificity and Matthew's 
correlation coefficient (MCC): 

tp 

sens = — (14) 



tp + fn 
tn 

tn + fp 



tn 

spec = — (15) 



tp-tn-fp-fn 
^{tp + fp) ■ {tp + fn) ■ {tn + fp) ■ {tn + fn) 

Sensitivity and specificity are always between zero and one and relate to how well the algorithms are 
able to find genes (true positives) and non-coding regions (true negatives), respectively. MCC reflects 
the overall correctness and lies between —1 and 1, where 1 represents perfect prediction. 

Gene level. When looking at the decoding position by position, genes that are predicted correctly do 
not contribute to the measures in an equal manner, but rather, the longer the gene, the more contribution 
it brings. However, it is interesting how well the genes are recovered, independent of how long they are. 
To measure this, we consider a predicted gene as one true positive if it overlaps by at least 50% of a 
true gene and as one false positive if there is no true gene with which it overlaps by at least 50%. Each 
true gene for which there is no predicted gene that overlaps by at least 50% counts as one false negative; 
see Figure 4. In this context, true negatives, i.e., the areas where there was no gene and no gene was 
predicted, are not considered, as they are not informative. The final measures are the recall, precision 
and F- score: 
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rec ■ prec 

f-score = 2 (19) 

rec + prec 

(20) 



2 ■ tp + fn + fp 



They are all between zero and one and reflect how well the true genes have been recovered. The 
recall gives the fraction of true genes that have been found, while the precision gives the fraction of the 
predicted genes that are true genes. The F-score is the harmonic mean of the two. 

Figure 4. Error types at the gene level. A predicted gene is considered one true positive if 
it overlaps with at least 50% of a true gene and one false positive if there is no true gene 
with which it overlaps by at least 50%. Each true gene for which there is no predicted gene 
that overlaps by at least 50% counts as one false negative. True genes that are covered more 
than 50% by predicted genes, but for which there is no single predicted gene that covers a 
minimum of 50% are disregarded. 



true — [I 



decoding — | | 1 || 1 

tp tp fn fp fn fp 



Figures 5 and 6 show the quality of predictions at the nucleotide and gene level, respectively. When 
comparing the Viterbi algorithm with the restricted Viterbi algorithm, it is clear that the restriction 
brings a great improvement to the prediction, as the restricted version has an increased power in all 
measures considered, with the exception of precision, where the Viterbi algorithm shows a tendency of 
increased precision with sequence length. However, when comparing the posterior- Viterbi algorithm 
with its restricted version, it is apparent that the restricted version does as good at the nucleotide level, 
but it performs worse at the gene level. By inspecting the predictions of the two methods, it was clear 
that the restricted posterior- Viterbi obtained an increased number of genes dictated by the expectation 
by just fractioning the genes predicted by the posterior- Viterbi. We believe this happens because the 
posterior- Viterbi algorithm finds the best local decoding, and therefore, adding global information, such 
as the total number of genes in the prediction, does not aid in the decoding. On the other hand, as 
the Viterbi algorithm finds the best global decoding, using the extra information results in a significant 
improvement of the prediction. Overall, at the nucleotide level, the posterior- Viterbi shows the best 
performance, while at the gene level, the restricted Viterbi has the highest quality. One might expect 
this, given the nature of the algorithms. 

Apart from these experiments, we also ran the algorithms on the annotated E. coli genome (GenBank 
accession number U00096). In this set of experiments, we split the genome into sequences of a length 
of approximately 10,000, for which we ran the algorithms as previously described. We found the same 
trends as in the performance for simulated data (results not shown). When using the E. coli data, we 
also recorded the running time of the algorithms, and we found that the restricted algorithms are about 
45 times slower than the standard algorithms, which is faster than the worst case scenario, which would 
lead to a A; ■ \Q\'^ = 7-5^ = 175 slowdown, as the average expectation of the number of patterns per 
sequence was k = 7. 
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Figure 5. Prediction quality at the nucleotide level given by average sensitivity, specificity 
and Matthew's correlation coefficient (MCC) for the decoding algorithms. We ran the 
restricted decoding algorithms using the expectation calculated from the distribution returned 
by the restricted forward algorithm. The plots show a zoom-in of the three measures. As both 
sensitivity and specificity are between zero and one, the Y-axes in these two plots have the 
same scale. 
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Figure 6. Prediction quality at the gene level given by average recall, precision and F-score 
for the decoding algorithms. We ran the restricted decoding algorithms using the expectation 
calculated from the distribution returned by the restricted forward algorithm. The plots show 
a zoom-in of the three measures with the same scale on the Y axes. 
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4. Conclusions 



We have introduced three novel algorithms that efficiently combine the theory of Hidden Markov 
Models with automata and pattern matching to recover pattern occurrences in the hidden sequence. 
First, we computed the distribution of the number of pattern occurrences by using an algorithm similar 
to the forward algorithm. This problem has been treated in [16] by means of Markov chain embedding, 
using simple finite sets of strings as patterns. Our method is, however, more general, as it allows the use 
of regular expressions. 
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From the occurrence number distribution, we calculated the expected number of pattern matches, 
which estimated the true number of occurrences with high precision. We then used the distribution 
to alter the prediction given by the two most widely used decoding algorithms: the Viterbi algorithm 
and the posterior- Viterbi algorithm. We have shown that in the case of the Viterbi algorithm, which 
finds the best global prediction, using the expected number of pattern occurrences greatly improves the 
prediction, both at the nucleotide and gene level. However, in the case of the posterior- Viterbi algorithm, 
which finds the best local prediction, such an addition only fragments the predicted genes, leading to a 
poorer prediction. Overall, deciding which algorithm is best depends on the final measure used, but as 
our focus was on finding genes, we conclude that the restricted Viterbi algorithm showed the best result. 

As the distribution obtained from the restricted forward algorithm facilitates the calculation of the 
distribution of the waiting time until the occurrence of the /c"" pattern match, the restricted Viterbi 
algorithm could potentially be further extended to incorporate this distribution while calculating the 
joint probability of observed and hidden sequences. Weighted transducers [26] are sequence modeling 
tools similar to HMMs, and analyzing patterns in the hidden sequence can potentially also be done by 
composition of the transducers, which describe the HMM and the automaton. 

Our method can presumably be used with already existing HMMs to improve their prediction, by 
using patterns that reflect the problems studied using the HMMs. For example, in [13], an HMM is 
used for finding frameshifts in coding regions. In this situation, the pattern would capture the sequence 
of hidden states that corresponds to a frameshift. In HMMs used for phylogenetic analysis [11,12], 
the hidden states represent trees, and an event of interest is a shift in the tree. A pattern capturing a 
change in the hidden state could thus aid in improving the prediction. In HMMs built for probabilistic 
alignments [10], the pattern could capture the start of an indel, and our method could potentially aid in 
finding a more accurate indel rate estimate. There is therefore a substantial potential in the presented 
method to successfully improve the power of HMMs. 
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